# Run regression
ipl12Reg <- lm(ipl12s ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
ipl24Reg <- lm(ipl24s ~ survey + immiFactor + gender+ yschoolingStd + ageStd + yearsInStd + sharedL, data = df)


ipl12SE <- coeftest(ipl12Reg, vcov= vcovHC(ipl12Reg, type="HC1"))
ipl24SE <- coeftest(ipl24Reg, vcov= vcovHC(ipl24Reg, type="HC1"))

# data frame with estimates
dfIpl12Reg <- data.frame(Variable = rownames(summary(ipl12Reg)$coef),
                         Coefficient = summary(ipl12Reg)$coef[, 1],
                         SE = ipl12SE[, 2],
                         modelName = "IPL-12")

dfIpl24Reg <- data.frame(Variable = rownames(summary(ipl24Reg)$coef),
                         Coefficient = summary(ipl24Reg)$coef[, 1],
                         SE = ipl24SE[, 2],
                         modelName = "IPL-24")

# levels(dfIpl12Reg$Variable)
# levels(dfIpl24Reg$Variable)

dfIpl12Reg$Variable <- factor(dfIpl12Reg$Variable,levels(dfIpl12Reg$Variable)[c(1, 7, 9, 8, 4, 5, 10, 6, 11, 3, 2)])
dfIpl24Reg$Variable <- factor(dfIpl24Reg$Variable,levels(dfIpl24Reg$Variable)[c(1, 7, 9, 8, 4, 5, 10, 6, 11, 3, 2)])

stargazer(ipl12Reg, ipl24Reg,
          se = list(ipl12SE[, 2],
                    ipl24SE[, 2]),
          digits = 3,
          align = T,
          no.space = TRUE,
          title = c("Full regression output for coefficent plots"),
          label = c("IPLRegs"), 
          covariate.labels = c("Sample B (Germany)", "Sample C (New York)", 
                               "Sample D (San Jose)", "Permanent Resident", 
                               "Naturalized", "Female", "Education (+1SD; 5y)",
                               "Age (+1SD; 18y)", "Residency (+1SD; 19y)", 
                               "Shared Language", "Constant"),
          dep.var.labels = c("IPL-12 Score rescaled", "IPL-24 Score rescaled"),
          notes = c("Robust standard errors in parentheses", 
                    "Reference sample: sample A (U.S. high income)", 
                    "Reference immigration status: temporary status"),
          out = c(here::here("Draft/PNAS/draft/tables", "IPLRegs.tex"))
)





# Further regressions (not in paper) after editor comments -----------------

# econ12sReg <- lm(econ12s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# econ12sRegSE <- coeftest(econ12sReg, vcov= vcovHC(econ12sReg, type="HC1"))
# 
# soc12sReg <- lm(soc12s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# soc12sRegSE <- coeftest(soc12sReg, vcov= vcovHC(soc12sReg, type="HC1"))
# 
# psy12sReg <- lm(psy12s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# psy12sRegSE <- coeftest(psy12sReg, vcov= vcovHC(psy12sReg, type="HC1"))
# 
# ling12sReg <- lm(ling12s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# ling12sRegSE <- coeftest(ling12sReg, vcov= vcovHC(ling12sReg, type="HC1"))
# 
# nav12sReg <- lm(nav12s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# nav12sRegSE <- coeftest(nav12sReg, vcov= vcovHC(nav12sReg, type="HC1"))
# 
# pol12sReg <- lm(pol12s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# pol12sRegSE <- coeftest(pol12sReg, vcov= vcovHC(pol12sReg, type="HC1"))
# 
# 
# stargazer(econ12sReg, soc12sReg, psy12sReg, ling12sReg, nav12sReg, pol12sReg,
#           se = list(econ12sRegSE[, 2],
#                     soc12sRegSE[, 2],
#                     psy12sRegSE[, 2],
#                     ling12sRegSE[, 2],
#                     nav12sRegSE[, 2],
#                     pol12sRegSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("IPL-12 Regs by Dimension"),
#           label = c("IPL12dims"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "IPL12dims.tex"))
# )
# 
# 
# 
# econ24sReg <- lm(econ24s ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# econ24sRegSE <- coeftest(econ24sReg, vcov= vcovHC(econ24sReg, type="HC1"))
# 
# soc24sReg <- lm(soc24s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# soc24sRegSE <- coeftest(soc24sReg, vcov= vcovHC(soc24sReg, type="HC1"))
# 
# psy24sReg <- lm(psy24s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# psy24sRegSE <- coeftest(psy24sReg, vcov= vcovHC(psy24sReg, type="HC1"))
# 
# ling24sReg <- lm(ling24s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# ling24sRegSE <- coeftest(ling24sReg, vcov= vcovHC(ling24sReg, type="HC1"))
# 
# nav24sReg <- lm(nav24s ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# nav24sRegSE <- coeftest(nav24sReg, vcov= vcovHC(nav24sReg, type="HC1"))
# 
# pol24sReg <- lm(pol24s  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# pol24sRegSE <- coeftest(pol24sReg, vcov= vcovHC(pol24sReg, type="HC1"))
# 
# 
# stargazer(econ24sReg, soc24sReg, psy24sReg, ling24sReg, nav24sReg, pol24sReg,
#           se = list(econ24sRegSE[, 2],
#                     soc24sRegSE[, 2],
#                     psy24sRegSE[, 2],
#                     ling24sRegSE[, 2],
#                     nav24sRegSE[, 2],
#                     pol24sRegSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("IPL-24 Regs by Dimension"),
#           label = c("IPL24dims"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "IPL24dims.tex"))
# )
# 
# # POL
# polUnderstanding_reg <- lm(polUnderstanding  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# polUnderstanding_regSE <- coeftest(polUnderstanding_reg, vcov= vcovHC(polUnderstanding_reg, type="HC1"))
# 
# polDiscuss_reg <- lm(polDiscuss  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# polDiscuss_regSE <- coeftest(polDiscuss_reg, vcov= vcovHC(polDiscuss_reg, type="HC1"))
# 
# polQuiz_reg <- lm(polQuiz  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# polQuiz_regSE <- coeftest(polQuiz_reg, vcov= vcovHC(polQuiz_reg, type="HC1"))
# 
# polAction_reg <- lm(polAction  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# polAction_regSE <- coeftest(polAction_reg, vcov= vcovHC(polAction_reg, type="HC1"))
# 
# 
# stargazer(polUnderstanding_reg, polDiscuss_reg, polQuiz_reg, polAction_reg,
#           se = list(polUnderstanding_regSE[, 2],
#                     polDiscuss_regSE[, 2],
#                     polQuiz_regSE[, 2],
#                     polAction_regSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("PolInt by Question"),
#           label = c("polInt"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "polInt.tex"))
# )
# 
# # ECON
# income_reg <- lm(income  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# income_regSE <- coeftest(income_reg, vcov= vcovHC(income_reg, type="HC1"))
# 
# employment3cat_reg <- lm(employment3cat  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# employment3cat_regSE <- coeftest(employment3cat_reg, vcov= vcovHC(employment3cat_reg, type="HC1"))
# 
# expenses_reg <- lm(expenses  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# expenses_regSE <- coeftest(expenses_reg, vcov= vcovHC(expenses_reg, type="HC1"))
# 
# emplSatisfied_reg <- lm(emplSatisfied  ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# emplSatisfied_regSE <- coeftest(emplSatisfied_reg, vcov= vcovHC(emplSatisfied_reg, type="HC1"))
# 
# stargazer(income_reg, employment3cat_reg, expenses_reg, emplSatisfied_reg,
#           se = list(income_regSE[, 2],
#                     employment3cat_regSE[, 2],
#                     expenses_regSE[, 2],
#                     emplSatisfied_regSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("EconInt by Question"),
#           label = c("econInt"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "econInt.tex"))
# )
# 
# # Soc
# dinner_reg <- lm(dinner ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# dinner_regSE <- coeftest(dinner_reg, vcov= vcovHC(dinner_reg, type="HC1"))
# 
# contactsInUsCat2_reg <- lm(contactsInUsCat2 ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# contactsInUsCat2_regSE <- coeftest(contactsInUsCat2_reg, vcov= vcovHC(contactsInUsCat2_reg, type="HC1"))
# 
# maxOrg_reg <- lm(maxOrg ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# maxOrg_regSE <- coeftest(maxOrg_reg, vcov= vcovHC(maxOrg_reg, type="HC1"))
# 
# favors_reg <- lm(favors ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# favors_regSE <- coeftest(favors_reg, vcov= vcovHC(favors_reg, type="HC1"))
# 
# stargazer(dinner_reg, contactsInUsCat2_reg, maxOrg_reg, favors_reg,
#           se = list(dinner_regSE[, 2],
#                     contactsInUsCat2_regSE[, 2],
#                     maxOrg_regSE[, 2],
#                     favors_regSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("SocInt by Question"),
#           label = c("socnInt"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "socInt.tex"))
# )
# 
# # Psy
# connectionHc_reg <- lm(connectionHc ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# connectionHc_regSE <- coeftest(connectionHc_reg, vcov= vcovHC(connectionHc_reg, type="HC1"))
# 
# outsider_reg <- lm(outsider ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# outsider_regSE <- coeftest(outsider_reg, vcov= vcovHC(outsider_reg, type="HC1"))
# 
# futureHc_reg <- lm(futureHc ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# futureHc_regSE <- coeftest(futureHc_reg, vcov= vcovHC(futureHc_reg, type="HC1"))
# 
# isolation_reg <- lm(isolation ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# isolation_regSE <- coeftest(isolation_reg, vcov= vcovHC(isolation_reg, type="HC1"))
# 
# stargazer(connectionHc_reg, outsider_reg, futureHc_reg, isolation_reg,
#           se = list(connectionHc_regSE[, 2],
#                     outsider_regSE[, 2],
#                     futureHc_regSE[, 2],
#                     isolation_regSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("PsyInt by Question"),
#           label = c("psyInt"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "psyInt.tex"))
# )
# 
# # Ling
# read_reg <- lm(read ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# read_regSE <- coeftest(read_reg, vcov= vcovHC(read_reg, type="HC1"))
# 
# speak_reg <- lm(speak ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# speak_regSE <- coeftest(speak_reg, vcov= vcovHC(speak_reg, type="HC1"))
# 
# write_reg <- lm(write ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# write_regSE <- coeftest(write_reg, vcov= vcovHC(write_reg, type="HC1"))
# 
# listen_reg <- lm(listen ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# listen_regSE <- coeftest(listen_reg, vcov= vcovHC(listen_reg, type="HC1"))
# 
# stargazer(read_reg, speak_reg, write_reg, listen_reg,
#           se = list(read_regSE[, 2],
#                     speak_regSE[, 2],
#                     write_regSE[, 2],
#                     listen_regSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("LingInt by Question"),
#           label = c("lingInt"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "lingInt.tex"))
# )
# 
# 
# # Nav
# seeDoctor_reg <- lm(seeDoctor ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# seeDoctor_regSE <- coeftest(seeDoctor_reg, vcov= vcovHC(seeDoctor_reg, type="HC1"))
# 
# findJobs_reg <- lm(findJobs ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# findJobs_regSE <- coeftest(findJobs_reg, vcov= vcovHC(findJobs_reg, type="HC1"))
# 
# legalProblems_reg <- lm(legalProblems ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# legalProblems_regSE <- coeftest(legalProblems_reg, vcov= vcovHC(legalProblems_reg, type="HC1"))
# 
# AdDTM_reg <- lm(AdDTM ~ survey + immiFactor + gender+ yschoolingStd + ageStd+ yearsInStd + sharedL, data = df)
# AdDTM_regSE <- coeftest(AdDTM_reg, vcov= vcovHC(AdDTM_reg, type="HC1"))
# 
# stargazer(seeDoctor_reg, findJobs_reg, legalProblems_reg, AdDTM_reg,
#           se = list(seeDoctor_regSE[, 2],
#                     findJobs_regSE[, 2],
#                     legalProblems_regSE[, 2],
#                     AdDTM_regSE[, 2]),
#           digits = 3,
#           align = T,
#           no.space = TRUE,
#           title = c("NavInt by Question"),
#           label = c("navInt"), 
#           font.size = "scriptsize",
#           keep.stat = c('n', 'adj.rsq'),
#           out = c(here::here("Draft/PNAS/draft/tables", "navInt.tex"))
# )






